/*
 * main.cpp
 *
 *  Created on: June 7, 2011
 *      Author: gibson
 */

#include <time.h>
#include <cmath>
#include <iostream>
#include "DataReader.h"
#include "TCanvas.h"
#include "TEnv.h"
#include "TApplication.h"
#include "TFile.h"
#include "TH2D.h"
#include "TH1D.h"
#include "TColor.h"
#include "TFitResultPtr.h"
#include "TFitResult.h"
#include "TF1.h"
#include "globals.hh"

using namespace std;

const int numCrystals = 2716;

int numCanvas = 0;

void TCanvas::Init() {
	// Initialize the TCanvas members. Called by all constructors.

	// Make sure the application environment exists. It is need for graphics
	// (colors are initialized in the TApplication ctor).
	if (!gApplication)
		TApplication::CreateApplication();

	// Get some default from .rootrc. Used in fCanvasImp->InitWindow().
	fHighLightColor = gEnv->GetValue("Canvas.HighLightColor", kRed);
	SetBit(kMoveOpaque, gEnv->GetValue("Canvas.MoveOpaque", 0));
	SetBit(kResizeOpaque, gEnv->GetValue("Canvas.ResizeOpaque", 0));
	if (gEnv->GetValue("Canvas.ShowEventStatus", kFALSE))
		SetBit(kShowEventStatus);
	if (gEnv->GetValue("Canvas.ShowToolTips", kFALSE))
		SetBit(kShowToolTips);
	if (gEnv->GetValue("Canvas.ShowToolBar", kFALSE))
		SetBit(kShowToolBar);
	if (gEnv->GetValue("Canvas.ShowEditor", kFALSE))
		SetBit(kShowEditor);
	if (gEnv->GetValue("Canvas.AutoExec", kTRUE))
		SetBit(kAutoExec);

	// Fill canvas ROOT data structure
	fXsizeUser = 0;
	fYsizeUser = 0;
	fXsizeReal = 20;
	fYsizeReal = 20;

	fDISPLAY = "$DISPLAY";
	fUpdating = kFALSE;
	fRetained = kTRUE;
	fSelected = 0;
	fClickSelected = 0;
	fSelectedX = 0;
	fSelectedY = 0;
	fSelectedPad = 0;
	fClickSelectedPad = 0;
	fPadSave = 0;
	fEvent = -1;
	fEventX = -1;
	fEventY = -1;
	fContextMenu = 0;

	numCanvas++;
}

void TCanvas::Closed() {
	Emit("Closed()");

	numCanvas--;
	if (numCanvas == 0)
		gApplication->Terminate(0);
}

double getCorrCoef(int* streamSigned, unsigned int numSamples, int pedestal) {
	unsigned short min = 1024;
	unsigned short max = 0;
	for (int i = 0; i < numSamples; i++) {
		if (streamSigned[i] < min)
			min = streamSigned[i];
		if (streamSigned[i] > max)
			max = streamSigned[i];
	}
	if ((max - min) > 64)
		return -2;

	if (pedestal != 0)
		pedestal = min;

	unsigned short* stream = new unsigned short[numSamples];
	for (int i = 0; i < numSamples; i++) {
		//		stream[i] = (unsigned short) streamSigned[i];
		if (streamSigned[i] > pedestal)
			stream[i] = (unsigned short) streamSigned[i] - (unsigned short) pedestal;
		else
			stream[i] = 0;
	}

	char const referenceSize = 21;
	unsigned long refTimesSize[referenceSize] = { 27, 77, 196, 381, 615, 854, 1056, 1188, 1230, 1180, 1053, 859, 659,
			476, 325, 211, 132, 86, 61, 44, 42 };

	unsigned long dataSum = 0;
	unsigned long data2Sum = 0;
	unsigned long dataRefSum;

	for (int i = 0; i < referenceSize; i++) {
		dataSum += (stream)[i];
		data2Sum += (stream)[i] * (stream)[i];
	}

	unsigned long const NSMRS = numSamples - referenceSize;
	unsigned long const RSTR2SMOTT = 364;

	double coefficient = 0;
	for (int i = 0; i < NSMRS; i++) {
		unsigned short* const tempStream = &(stream)[i];

		unsigned long DRSTRS = tempStream[0] * refTimesSize[0] + tempStream[1] * refTimesSize[1] + tempStream[2]
				* refTimesSize[2] + tempStream[3] * refTimesSize[3] + tempStream[4] * refTimesSize[4] + tempStream[5]
				* refTimesSize[5] + tempStream[6] * refTimesSize[6] + tempStream[7] * refTimesSize[7] + tempStream[8]
				* refTimesSize[8] + tempStream[9] * refTimesSize[9] + tempStream[10] * refTimesSize[10]
				+ tempStream[11] * refTimesSize[11] + tempStream[12] * refTimesSize[12] + tempStream[13]
				* refTimesSize[13] + tempStream[14] * refTimesSize[14] + tempStream[15] * refTimesSize[15]
				+ tempStream[16] * refTimesSize[16] + tempStream[17] * refTimesSize[17] + tempStream[18]
				* refTimesSize[18] + tempStream[19] * refTimesSize[19] + tempStream[20] * refTimesSize[20];

		if (DRSTRS > (dataSum << 9)) {
			unsigned long DRSTRSMDS = DRSTRS - (dataSum << 9);
			double coefTemp = sqrt(
					(double) (DRSTRSMDS * DRSTRSMDS) / (double) ((RSTR2SMOTT * (referenceSize * data2Sum - dataSum
							* dataSum))) / 512);
			if (coefTemp > coefficient) {
				coefficient = coefTemp;
			}
		}

		dataSum += tempStream[referenceSize] - tempStream[0];
		data2Sum += tempStream[referenceSize] * tempStream[referenceSize] - tempStream[0] * tempStream[0];
	}

	delete[] stream;

	return coefficient;
}

int main(int argc, char** argv) {
	int argcFake = 0;
	TApplication* app = new TApplication("app", &argcFake, argv);

	string baseDir = "";

	int runNum = 0;
	bool expand = true;
	bool view = false;
	bool write = true;
	string output = "./data/histo.root";
	for (int i = 1; i < argc; i++) {
		if (argv[i][0] == '-') {
			if (argv[i][1] == 'd') {
				baseDir = argv[i + 1];
			} else if (argv[i][1] == 'r')
				runNum = atoi(argv[i + 1]);
			else if (argv[i][1] == 's')
				expand = false;
			else if (argv[i][1] == 'V') {
				Double_t r[5] = { 1, 1, 0, 0, 0 };
				Double_t g[5] = { 0, 1, 1, 1, 0 };
				Double_t b[5] = { 0, 0, 0, 1, 1 };
				Double_t stop[5] = { 0, .25, .5, .75, 1 };
				TColor::CreateGradientColorTable(5, stop, r, g, b, 20);

				TFile* f = new TFile(argv[i + 1]);
				TH2D* he = (TH2D*) f->Get("histoE");
				TH2D* hp = (TH2D*) f->Get("histoP");
				TH1D* hz = (TH1D*) f->Get("histoZ");
				TH1D* he1 = (TH1D*) f->Get("histoE1");
				TH1D* hc = (TH1D*) f->Get("histoC");
				TH1D* hl = (TH1D*) f->Get("histoL");
				TH2D* hph6 = (TH2D*) f->Get("histoPh6");
				TH2D* hph65 = (TH2D*) f->Get("histoPh65");
				TH2D* hph7 = (TH2D*) f->Get("histoPh7");
				TH1D* hphE = (TH1D*) f->Get("histoPhEne");
				TCanvas* c1 = new TCanvas("c1");
				he->Draw("colz");
				TCanvas* c2 = new TCanvas("c2");
				hp->Draw("colz");
				TCanvas* c3 = new TCanvas("c3");
				he1->Draw();
				TCanvas* c4 = new TCanvas("c4");
				hz->Draw();
				TCanvas* c5 = new TCanvas("c5");
				hc->Draw();
				TCanvas* c6 = new TCanvas("c6");
				hl->Draw();
				TCanvas* c7 = new TCanvas("c7", "c7", 800, 800);
				hph7->Draw("colz");
				TCanvas* c8 = new TCanvas("c8", "c8", 800, 800);
				hph65->Draw("colz");
				TCanvas* c9 = new TCanvas("c9", "c9", 800, 800);
				hph6->Draw("colz");
				TCanvas* c10 = new TCanvas("c10");
				hphE->Draw();
				app->Run(true);
				f->Close();
				app->Delete();

				return 0;
			} else if (argv[i][1] == 'v') {
				view = true;
			} else if (argv[i][1] == 'n' && argv[i][2] == 'w')
				write = false;
			else if (argv[i][1] == 'o')
				output = argv[i + 1];
		}
	}

	if (!baseDir.compare("")) {
		cout << "input base directory with -d flag" << endl;
		return -1;
	}

	DataReader reader(baseDir, runNum, expand);

	TFile* f = new TFile(output.c_str(), "recreate");
//	TH2D* he = new TH2D("2DEneF", "Corr vs Energy for non-zero energy channels;Correlation Coefficient;Crystal Energy(MeV)", 100, 0, 1, 200, 0, 4);
	TH2D* he = new TH2D("2DEneF", "Corr vs Energy for non-zero energy channels;Correlation Coefficient;Crystal Energy(MeV)", 100, 0, 1, 200, 0, 1);
	TH2D* he2 = new TH2D("2DEneF2", "Corr vs Energy including zero energy channels;Correlation Coefficient;Crystal Energy(MeV)",100,0,1,200,0,4);
	TH2D* hp = new TH2D("2DPed", "Correlation coefficient vs unaccounted for pedestal;Correlation Coefficient;Pedestal(count)", 100, 0, 1, 256, 0, 1023);
	TH1D* hz = new TH1D("1DZer", "Correlation coefficient for 0 energy crystals;Correlation Coefficient;Channel Count", 100, 0, 1);
	TH1D* he1 = new TH1D("1DEne", "Correlation coefficient for 1 MeV crystals;Correlation Coefficient;Channel Count", 100, 0, 1);
	TH1D* hc = new TH1D("1DCom", "% Compression by corr;Correlation Coefficient;Compression(%)", 100, 0, 1);
	TH1D* hl = new TH1D("1DLos", "% of a photon's energy passed by discrimination by corr;Correlation Coefficient; Saved Energy(%)", 100, 0, 1);
//	TH2D* hph6 = new TH2D("2DPhot6", "% of a photon's energy passed by discrimination with a 0.6 threshold;Photon Energy(MeV);Saved Energy(%)", 200, 0, 200, 200, 0, 100);
//	TH2D* hph65 = new TH2D("2DPhot65", "% of a photon's energy passed by discrimination with a 0.65 threshold;Photon Energy(MeV);Saved Energy(%)", 200, 0, 200, 200, 0, 100);
//	TH2D* hph7 = new TH2D("2DPhot7", "% of a photon's energy passed by discrimination with a 0.7 threshold;Photon Energy(MeV);Saved Energy(%)", 200, 0, 200, 200, 0, 100);
//	TH1D* hphE = new TH1D("1DPhotEne", "Photon Energy Distribution;Photon Energy(MeV);Count", 1500, 0, 1500);
	TH2D* hph6 = new TH2D("2DPhot6", "% of a photon's energy passed by discrimination with a 0.6 threshold;Photon Energy(MeV);Saved Energy(%)", 55, .5, 5, 200, 0, 100);
	TH2D* hph65 = new TH2D("2DPhot65", "% of a photon's energy passed by discrimination with a 0.65 threshold;Photon Energy(MeV);Saved Energy(%)", 50, 0.5, 1, 200, 0, 100);
	TH2D* hph7 = new TH2D("2DPhot7", "% of a photon's energy passed by discrimination with a 0.7 threshold;Photon Energy(MeV);Saved Energy(%)", 50, .5, 1, 200, 0, 100);
	TH1D* hphE = new TH1D("1DPhotEne", "Photon Energy Distribution;Photon Energy(MeV);Count", 50, .5, 1);

	double energyTot = 0;
	double energy[2760];
	int** dataTemp = reader.GetData(0, &energy[0]);
	int l;
	for (l = 0; true; l++) {
		cout << "\r" << l << flush;
		double energyPhot = 0;
		double savedEnergy6 = 0;
		double savedEnergy65 = 0;
		double savedEnergy7 = 0;
		for (int j = 0; j < numCrystals; j++) {
			if (dataTemp[j] != NULL) {
				energyPhot += energy[j];
				double coef = getCorrCoef(dataTemp[j], 48, dataTemp[j][50]);
				if (coef == -2) {
					savedEnergy7 += energy[j];
					savedEnergy65 += energy[j];
					savedEnergy6 += energy[j];

					for (double k = .005; k <= 1; k += .01){
						hc->Fill(k);
						hl->Fill(k,energy[j]);
					}

					delete[] dataTemp[j];
					continue;
				}

				if (coef >= 0.7) {
					savedEnergy7 += energy[j];
				}
				if (coef >= 0.65) {
					savedEnergy65 += energy[j];
				}
				if (coef >= 0.6) {
					savedEnergy6 += energy[j];
				}

				if (energy[j] > 0) {
					hp->Fill(getCorrCoef(dataTemp[j], 48, 0), dataTemp[j][50]);
				} else if (energy[j] == 0)
					hz->Fill(coef);
				if (energy[j] > .9 && energy[j] < 1.1) {
					he1->Fill(coef);
				}
				if(energy[j]>0)
					he->Fill(coef,energy[j]);
				he2->Fill(coef, energy[j]);
				for (double k = .005; k <= coef; k += .01) {
					hc->Fill(k);
					hl->Fill(k, energy[j]);
				}
			}

			delete[] dataTemp[j];
		}
		delete[] dataTemp;
		dataTemp = reader.GetData(l + 1, &energy[0]);
		if (dataTemp == NULL) {
			break;
		}

		hph7->Fill(energyPhot, savedEnergy7 / energyPhot * 100);
		hph65->Fill(energyPhot, savedEnergy65 / energyPhot * 100);
		hph6->Fill(energyPhot, savedEnergy6 / energyPhot * 100);

		hphE->Fill(energyPhot);

		energyTot += energyPhot;
	}
	cout << endl;

	for (int i = 1; i <= hc->GetNbinsX(); i++) {
		hc->SetBinContent(i, (1 - (hc->GetBinContent(i) / ((l+1) * numCrystals))) * 100);
		}
	for (int i = 1; i <= hl->GetNbinsX(); i++) {
		hl->SetBinContent(i, (hl->GetBinContent(i) / energyTot) * 100);
	}

	Double_t r[5] = { 1, 1, 0, 0, 0 };
	Double_t g[5] = { 0, 1, 1, 1, 0 };
	Double_t b[5] = { 0, 0, 0, 1, 1 };
	Double_t stop[5] = { 0, .25, .5, .75, 1 };
	TColor::CreateGradientColorTable(5, stop, r, g, b, 20);

	if (write) {
		he->Write("histoE");
		he2->Write("histoE2");
		hp->Write("histoP");
		he1->Write("histoE1");
		hz->Write("histoZ");
		hc->Write("histoC");
		hl->Write("histoL");
		hph7->Write("histoPh7");
		hph65->Write("histoPh65");
		hph6->Write("histoPh6");
		hphE->Write("histoPhEne");
	}
	if (view) {
		TCanvas* c1 = new TCanvas("c1");
		he->Draw("colz");
		TCanvas* c2 = new TCanvas("c2");
		hp->Draw("colz");
		TCanvas* c3 = new TCanvas("c3");
		he1->Draw();
		TCanvas* c4 = new TCanvas("c4");
		hz->Draw();
		TCanvas* c5 = new TCanvas("c5");
		hc->Draw();
		TCanvas* c6 = new TCanvas("c6");
		hl->Draw();
		TCanvas* c7 = new TCanvas("c7", "c7", 800, 800);
		hph7->Draw("colz");
		TCanvas* c8 = new TCanvas("c8", "c8", 800, 800);
		hph65->Draw("colz");
		TCanvas* c9 = new TCanvas("c9", "c9", 800, 800);
		hph6->Draw("colz");
		TCanvas* c10 = new TCanvas("c10");
		hphE->Draw();
		app->Run(true);
	}

	f->Close();
	app->Delete();

	return 0;
}

